.ipynb .pdf

repository open issue

Contents

# # From XML data
#
# Example on how to load your data if it is currently in xml format
import json

import dolfin
from fenics_plotly import plot

import pulse

try:
    from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
    from dolfin import Mesh, Function, Constant, DirichletBC, interpolate

# Mesh
mesh = Mesh(pulse.utils.mpi_comm_world(), "data/mesh.xml")

# Marker functions
facet_function = dolfin.MeshFunction("size_t", mesh, "data/facet_function.xml")
marker_functions = pulse.MarkerFunctions(ffun=facet_function)

# Markers
with open("data/markers.json", "r") as f:
    markers = json.load(f)

# Fiber
fiber_element = dolfin.VectorElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=4,
    quad_scheme="default",
)
fiber_space = dolfin.FunctionSpace(mesh, fiber_element)
fiber = Function(fiber_space, "data/fiber.xml")

microstructure = pulse.Microstructure(f0=fiber)

# Create the geometry
geometry = pulse.HeartGeometry(
    mesh=mesh,
    markers=markers,
    marker_functions=marker_functions,
    microstructure=microstructure,
)

activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
activation.assign(Constant(0.0))
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
    activation=activation,
    parameters=matparams,
    f0=geometry.f0,
)

# LV Pressure
lvp = Constant(0.0)
lv_marker = markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]

# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [pulse.RobinBC(value=Constant(base_spring), marker=markers["BASE"][0])]

# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    bc = DirichletBC(V.sub(0), Constant(0.0), geometry.ffun, markers["BASE"][0])
    return bc

dirichlet_bc = [fix_basal_plane]
# You can also use a built in function for this
# from functools import partial
# dirichlet_bc = partial(pulse.mechanicsproblem.dirichlet_fix_base_directional,
#                        ffun=geometry.ffun,
#                        marker=geometry.markers["BASE"][0])

# Collect boundary conditions
bcs = pulse.BoundaryConditions(
    dirichlet=dirichlet_bc,
    neumann=neumann_bc,
    robin=robin_bc,
)

# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

problem.solve()

(2, True)

# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, 0.2))

2021-11-23 08:30:29,488 [782] INFO     pulse.iterate: Iterating to target control (f_26)...

2021-11-23 08:30:29,491 [782] INFO     pulse.iterate: Current control: f_26 = 0.000,0.000

2021-11-23 08:30:29,492 [782] INFO     pulse.iterate: Target: 1.000,0.200

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 78),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 252),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 389),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 450),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 511),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 572),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 641),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 762),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 823),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 892),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1005),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1110),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1275),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1328),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1405),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1519),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1589)],
 [(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 74),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 76)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 248),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 250)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 385),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 387)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 446),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 448)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 507),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 509)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 568),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 570)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 637),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 639)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 758),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 760)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 819),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 821)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 888),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 890)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1001),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1003)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1106),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1108)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1271),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1273)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1324),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1326)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1401),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1403)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1515),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1517)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1585),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1587))])

# Get the solution
u, p = problem.state.split(deepcopy=True)

volume = geometry.cavity_volume(u=u)
print(f"Cavity volume = {volume}")

Cavity volume = 1.1409709245419595

# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, show=False))
fig.show()

previous

<no title>

next

<no title>